Phase Dynamics in Intrinsic Josephson Junctions and its Electrodynamics 
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We present a theoretical description of the phase dynamics and its corresponding electrodynamics 
in a stack of inductively coupled intrinsic Josephson junctions of layered high-Tc superconductors in 
the absence of an external magnetic field. Depending on the spatial structure of the gauge invariant 
phase difference, the dynamic state is classified into: state with kink, state without kink, and state 
with solitons. It is revealed that in the state with phase kink, the plasma is coupled to the cavity 
and the plasma oscillation is enhanced. In contrast, in the state without kink, the plasma oscillation 
is weak. It points a way to enhance the radiation of electromagnetic from high- Tc superconductors. 
We also perform numerical simulations to check the theory and a good agreement is achieved. The 
radiation pattern of the state with and without kink is calculated, which may serve as a fingerprint 
of the dynamic state realized by the system. At last, the power radiation of the state with solitons 
is calculated by simulations. The possible state realized in the recent experiments is discussed in the 
viewpoint of the theoretical description. The state with kink is important for applications including 
terahertz generators and amplifiers. 

PACS numbers: 74.50.+r, 74.25.Gz, 85.25.Cp 



I. INTRODUCTION 

The electromagnetic waves in the terahertz(THz) re- 
gion, which is defined in the range from O.lTHz to lOTHz, 
have wide applications, such as drug detection, materi- 
als characterization, security check, and so on. This has 
sparked considerable efforts to seek compact and low-cost 
solid-state generators. [J, 

It has been known for a long time that Josephson junc- 
tions can be used as electromagnetic oscillators. The 
power radiated from a single junction, however is in 
the range of pW, which is too small for practical ap- 
plications. The frequency is about one hundred giga- 
hertz because of the small superconducting energy gap 
for conventional superconductors. [1, H, [1, 01 Although 
one may integrate a large array of Josephson junctions 
made of conventional superconductors on a chip to en- 
hance the radiation power, 0, H, H, the frequency is 
still below terahertz. The discovery of intrinsic Joseph- 
son junctions in layered high- Tc superconductors, such 
as Bi2Sr2CaCu208+5(BSCCO) provides a very nice can- 
didate for terahertz oscillator. [llj The advantages of in- 
trinsic Josephson junctions over conventional low tem- 
perature junctions are as follows. First, the junctions 
are homogeneous in the atomic scale, which makes the 
coherent radiation in large number of junctions possible. 
Secondly, the energy gap is about 60meV which corre- 
sponds to 15THz. The terahertz Josephson plasma if 
excited is thus free from Landau damping. [T^ 

One idea to excite the terahertz wave inside the in- 
trinsic Josephson junctions is by the motion of Joseph- 
son vortices lattice induced by an in-plane magnetic field 
and a transport current, which has already been in- 
vestigated both theoretically and experimentally. [H, [3, 
m, [M [13, H Ell In spite of these works, it is still 



lack of clear evidence of coherent radiation. Radiation 
from BSCCO with iiijection of quasiparticles has been 
reportedfH iU, [H |2l, [2l|. Alternatively, the tera- 
hertz radiation without a magnetic field has also been 
attempted m, [H, [13, [H, HI- Recently, a strong co- 
herent radiation from BSCCO in the absence of a mag- 
netic field was observed [2?! [29| . where the mesa of the 
single crystal of BSCCO forms a cavity. The break- 
through in the experiments has inspired considerable the- 
oretical and experimental efforts, aiming to reveal the 
mechanism of strong radiation, [s^, HU, [H, [H, [s^, [111 
A new dynamic state has been suggested to explain the 
experiments 0, |3l| . 

It is well known that there exist various dynamic 
states, such as the McCumber state and states with soli- 
tons, with different IV characteristics in a Josephson 
iunction[30:, ,36,] due to the nonlinearity. In the McCum- 
ber state, the gauge invariant phase difference is uniform 
in space, which we will refer to as the state without kink 
in later discussions. The soliton solutions are also well 
known, especially in a single junction, where a quantized 
particle-like object of 27r phase variation travels along 
the junction. Recently, a new dynamic state was found, 
where a (2m-|- l)7r phase kink is localized inside the junc- 
tion with an integer m. [IS [111] We will refer to this state 
hereafter as state with kink. Because of the complexity 
in the dynamics in this highly nonlinear system, theoret- 
ical understanding of phase dynamics and its electrody- 
namics, and finding the optimal state are expected to be 
helpful for experimental realizations of terahertz gener- 
ators. Assessments of such a device from a theoretical 
aspect are also needed. 

For this purpose, in this paper we provide more details 
on the new dynamic state found in the previous study. ^30j 
In Section II, we first derive the Lagrangian of a stack 
of Josephson junctions based on the superconductor- 
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FIG. 1: (Color online). Schematic view of a stack of Joseph- 
son junctions based on the SIS model. The blue(pink) area 
denotes superconducting(insulating) layers. 

insulator-superconductor(SIS) model. From the La- 
grangian we derive the inductively coupled sine-Gordon 
equation. It also gives the power balance condition. In 
section III, we develop a general procedure to solve the 
coupled sine-Gordon equation from the spectrum analy- 
sis, from which we derive the solutions with and without 
phase kink. We calculate analytically the IV character- 
istics and power radiation of the state with and without 
kink when the plasma oscillation is small. Numerical 
simulations are performed to verify the analytical solu- 
tions and a good agreement is attained. In section IV, 
the energy stored in the system is evaluated. In section 
V, the far-field radiation pattern from the mesa is calcu- 
lated both at state with and without kink, which can be 
used to distinguish the different dynamic states. In sec- 
tion VI, the power radiation from the state with solitons 
is investigated by simulations for comparison. At last, 
the paper is concluded with a short discussion. 

II. LAGRANGIAN AND MODEL EQUATION 

The geometry we consider is depicted in Fig. [TJ We 
neglect the thermal fluctuations so that the amplitude 



|A| of the superconducting order parameter |A| expiitp) 
is constant. Furthermore, ■0 along the y axis is assumed 
to be uniform, namely we concentrate on the zero mode 
along this direction. The system is then reduced to two 
dimensions. The density of energy stored in the super- 
conducting layers, which consists of supercurrent energy 
and magnetic energy, then can be written as 



Hsiix) = - / [A?(V X B^Jl + ByJ]dz, (1) 

Jl(s+D)+D 



where is the penetration depth, and i?^; is the mag- 
netic field in the Ith. superconducting layer. The energy 
stored in the insulating layers is the sum of energy of 
electromagnetic wave and Josephson energy 



fl(s+D)+D oy2 ^ 

Ji{s+D) Stt Stt 27rc 

(2) 

where [Ei) is the magnetic (electric) field in the ^th 
insulating layer (their variation along the z direction in 
the ^th insulating layer will be neglected in later treat- 
ment). Jc is the critical current density, ec dielectric con- 
stant along the z axis and c the light velocity in vacuum. 
Pi is the gauge invariant phase difference defined as 



2^ Ms+D)+D 

Pi{x) - iJi+i{x) - V/(x) - — / At{x)dz, (3) 

^0 Jl(s+D) 



where $o = hc/2e is the flux quantum and Af 
is the vector potential. The magnetic field inside 
the superconducting layer can be evaluated from the 
London equation B^i{z) = (sinh[(s — z)/Xs]Bf + 
sinh[z/As]-Bf_^^)/ sinh[s/As]. In high- Tc superconductors, 
the thickness of superconducting layer s and insulat- 
ing layer D is much smaller than A^, we have B^^ ~ 
[(s — z)Bf + zBf^-^]/s. Then the total energy density can 
be expressed as 



H{x) = J2[Hsi{x) + Hu{x)] = E{^[(2C + l)Bf - aSl^Bf + Bf Bl,)] + d'-^ + ^ J,(l - cosP,)}, (4) 

I 

where we have neglected the surface effect along the stack be noted that the inductive coupling is very strong in 
direction, which is valid for thick stacks of junctions. ^ = BSCCO, see Table HI 
Xl/sD is the strength of inductive coupling. It should 
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TABLE I; Conversion of quantities among dimensionless, Gaussian and SI units. Here Ac and Xab are the penetration depth; 
Ec is dielectric constant along the z axis; c is the light velocity in vacuum; eg is dielectric constant in vacuum; uip = c/Xc^/e^ is 
the Josephson plasma frequency. Jc = c^o/^n^ X^D is the critical current density. In the present paper, we use Ac ~ 200/^m, 
Aab = OAfim, Ec ~ 10, s — 0.3nm and D — 1.2nm, which are typical for BSCCO. Then A^ — \J sD/[s + D)'^Xab ~ 0.16/im. The 
two dimensionless parameters are then /3 = 0.02 and C = 7.1 x lO"*. Following Ref. we use slightly larger = 4/9 x 10 in 
the present paper, and the results are insensitive to ^ as far as it is large. The length of junction is L = 80/im. 
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To find the relation between the magnetic field Bf and 
Pi, we derive both sides of Eq. ([3]) with respect to x. 
With the London equation and Maxwell equation, we 
arrive at 



^0 

2ttD 



d^P = MBy, 



(5) 



where P is a column vector = [Pi, P2, Pn] with N 
being the number of junctions. The column vectors for 
other quantities are defined in the same way. M is the 
inductive coupling matrix defined as[37l [38| 



M = 
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(6) 

where the periodic boundary condition along the z axis 
is imposed. Using the ac Josephson relation 9tP = 
2eW-D/h and Eq. (O, we can rewrite the total energy 
density in a compact form 

H{x) = ia,P^M-i9,P + i^tP^StP + ^(1 - cosP/), 

(7) 

where the dimensionless quantities have been used, 
which, with the conversion among SI units and Gaus- 
sian units, are compiled in Table HI The first term at the 
right-hand side of Eq. Q represents the magnetic energy, 
the second term the electric energy and the last term the 
Josephson coupling. The Lagrangian corresponding to 
Eq. is 

^{x) = i5,P^M-i5,P + i5*P^5tP-^(l-cosP,). 

I 

(8) 

With the Euler-Lagrangian formula, we arrive at the cou- 
pled sine- Gordon equation 



where sinP = [sin Pi^sin P2, sin Pjv]'^- We consider 
the overlap geometry [37[ where the current is uniformly 
injected into the system. Taking the dissipation and ex- 
ternal current into account, we obtain the inductively 
coupled perturbed sine-Gordon equation 



dlV = M[sinP 



■/3c)tP-Jext], 



(10) 



where the first term at the right-hand side of Eq. (fTO|) is 
the Josephson current, the second term the displacement 
current, the third term the quasiparticles contribution 
and the last term the external current. Writing down 
the equation for P; in Eq. pn]) explicitly, we have 

dlPi - (1 - CA(2) ) [sin Pi + d^tPi + l3dtPi - Jcxt] , (11) 

where A^^^ is the finite difference operator defined as 
A^^)// = fi+i + /;_! — 2//. Besides the inductive cou- 
pling in Eq. (fTU|) . a capacitive coupling [sol lioj and a 
coupling originating from non-equilibrium effect s[4ll 
are also present in intrinsic Josephson junctions. These 
two couplings are weak in comparison to the inductive 
coupling and are neglected in the present work. 

The above calculations are based on the SIS model 
(superconductor-insulator-superconductor), which is a 
good model for artificially stacked Josephson junctions. 
However, for BSCCO, the superconducting layer is only 
of atomic thickness, and thus the quantity is not well 
defined. To find the relation between and the mea- 
surable penetration depth Aab, we need to resort to the 
Lawrence-Doniach model, which has already been dis- 
cussed extensively in literatures. The connection between 
A, and Xab is given by A, = ^ sD / {s + DfXab.^M, 

In the presence of dissipations and a power input, the 
energy oscillates with time according to 

dtH{x) = d^t'P'^M.-^d^V + {d^P"^ + sinP^)5fP. (12) 



a^p = M[sinP + 9t2p] 



(9) With the help of Eq. (fTO|) . we have for the steady state 
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/ / dtHdxdt^ {dtP'^M-^d^P)\l^dt+ / dtP^J^^tdxdt-f] / dtP^dtPdxdt ^ 0, (13) 
Jo Jo Jo Jo Jo Jo Jo 



where L is the length of junctions and T is the period of 
plasma oscillation. Rewriting Eq. (jl3p in a more trans- 
parent form, we have the power balance condition 



r iE'^B)\^dt + LTElJ, 
Jo 



(3 1 i E^'Edxdt = 0, 
Jo Jo 

(14) 

where the first term is the power gain (loss) at edges 
due to irradiation (radiation). The second term is the 
input power with Ejc the dc electric field, and the last 
term is the energy loss due to dissipations. It should be 
remarked that Eq. p4)) is general and should be valid 
at different states. As will be shown later, the power 
balance relation is useful when the plasma oscillation is 
strong and the linear expansion fails. From Eq. (|14p and 
E = Eac + Edc, we can see that when the oscillation in 
electric field Eac is small, the IV curve is almost ohmic 
Edc ~ Jext//9- To have strong radiation, the oscillation of 
electric field in the junctions should be large. Therefore 
the optimal state for the radiation is a state having non- 
linear IV characteristics where a large part of the input 
power can be pumped into plasma oscillation. The prob- 
lem then boils down to finding such a state with highly 
nonlinear IV characteristics. 

The relation between the oscillating magnetic field and 
electric field at the edges of junctions is given by the 
boundary condition. The boundary condition depends on 
many effects such as distribution of the order parameter 
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FIG. 2: (Color online), (a), (b): Schematic view of two sim- 
plest configurations of the static phase Pi(x). (c), (e), (g) 
and (i): (2m + l)7r phase kink for m = 0, 1, 2 and 3 respec- 
tively, (d), (f), (h) and (J); their corresponding un-quantized 
static vortices. 



near edges, the geometry of the sample and the dielectric 
materials attached to the sample. In Refs. [3, HE], the 
dynamic boundary condition is derived from the elec- 
tromagnetic wave equations inside the dielectrics. In 
the present paper, we use an effective impedance as the 
boundary condition 



Eac/Bac = Z= |Z(c^)| exp(i0(w)), 



(15) 



where \Z\ and d are parameters, and lo is the frequency. 
This boundary condition is general and any other bound- 
ary condition can be casted into this form. It has been 
pointed out that there is a significant impedance mis- 
match between the intrinsic Josephson junctions and out- 
side space because of the small ratio between the thick- 
ness of the stack and the penetration depth Ac.^] This 
means \Z\ » 1, which is similar to the single junction 
case.[46| 

The radiation power counted by Poynting vector at one 
edge with an effective impedance becomes 



ELEacrfi. 
(16) 

For ease of theoretical calculation, we consider the sit- 
uation that the radiation does not substantially change 
the plasma oscillation inside the junctions. In this case, 
Eac can be evaluated without radiation, i.e. with the 
simple boundary condition dxP — 0. This treatment is 
valid when the impedance mismatch is significant, which 
will be shown later to be rather accurate by numerical 
simulations. When the impedance \Z\ is small, however, 
one should take the effect of radiation into account in the 
calculation of the plasma dynamics self-consistently. 



III. SOLUTIONS 

In this section, we first construct a general proce- 
dure to solve the coupled sine-Gordon equation from the 
spectrum analysis. There exist longitudinal and trans- 
verse plasma modes in a stack of junctions. As the 
stack itself forms a cavity, the plasma component can 
be written as Pi{x) cos(fcjx) sin(g//(iV -|- 1))[43, Sll 
with kj ~ JTr/i, and q, j integers. There are N dif- 
ferent dispersion branches with characteristic velocity 
Cq = l/v^l + 2C[1 - cos(q7r/(iV -|- 1))]. When the stack 
is thick enough, the plasma oscillation uniform along the 
c axis becomes possible and its velocity is cq = 1. We 
will concentrate on this case in the present work since it 
supports the strong radiation. Without an external in- 
plane magnetic field, the solution including all frequency 
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harmonics subject to the boundary condition dxPi = 
can be expressed as(30| 

oo 

Pi{x, t) = ujt + Pf{x) + ^ Re[—iAj exp(ijwt)] cos{kjx), 

(17) 

where ut is the rotating part, Pf the static phase kink 
and the last term is the plasma oscillation including 
all harmonics, with Aj the oscillation amplitude. For 
simplicity, the first cavity mode along the x axis with 



ki = tt/L is considered and the small time dependence of 



Pf is neglected. [301 131| Putting the mth frequency com- 
ponent of the Josephson current sinP; as Sm exp{imujt), 
and expanding sine of sine with Bessel functions, we have 



^0 



for m > 1 , 
for m = 0, 



(18) 



where 



4-00 



{(}j=-oo} j=l 



Y[jg^{\Aj\cos{kjx)) 



(19) 



-l-CXD 



with Jq. the Bessel function of the first kind. ^ is the summation over the ensemble of g^'s, and Aj — 

{<?3=-00} 

l^jl exp(i(/)j). Substituting Eq. (jl7p into the coupled sine-Gordon Eq. (|lip and comparing each frequency component 
in sinP;, we have for the mth (m > 1) component, 



[ikl, - i(mw)2 - Pmuj]A.^ cos(fc„a:) = (1 - C^'^'^^)SL 

I 



(20) 



From Eq. ([20)) . A„i is given by 

m. 



(21) 



with 



F„ 



- / (l-CA(2))5^„cos(fc™x)dx. (22) 
^ Jo 

The functional Fm represents the coupling of the plasma 
to the cavity modes, which is the central quantity for 
the excitation of plasma. Other mechanism such as the 
modulation of the critical current is also possible [28l.l49|. 
although it is practically very hard to achieve a homoge- 
neous modulation along the z axis. In the present solu- 
tion, Pf is inherently one part of the solution. Since we 
have assumed that phase is uniform along the y direc- 
tion, or equivalently, we have considered the (1, 0) cavity 
mode, only the kink in the x direction contributes to Fm- 
If there exits a kink along the y direction simultaneously, 
the functional F„i will be enhanced further, [s^ In the 
case of cylinder geometry, the plasma is coupled to the 
cavity fully via the kink so that F^ is maximized, fs^ 

It should be noted that Am is independent of / in Eq. 
(I2ip . which imposes a constraint on the arrangement of 
Pf in the z direction. As will be shown later, periodic 
arrangements such as those in Figs. Ufa) and (b) diag- 
onalize the finite difference operator and make Fm inde- 
pendent of 

For the 0th component (static part), we obtain 



(23) 



The current conservation relation reads Joxt = P'^+{Sq)x 
{{■■■)x is the spatial average). The remaining terms in Eq. 
which do not contribute to the net current is 



(24) 



where C >> 1 is taken into account. 

From Eqs. dTT]), ^ and we can calculate the 

plasma oscillation, IV characteristics and Pf. We con- 
sider explicitly the case where the fundamental mode 
Ai is small and thus higher harmonics can be safely 
neglected. We also approximate Ji(|Ai| cos(A;ia;)) ~ 
|Ai| cos(fcia;)/2 and Jo(| Ai | cos(fcix)) ~ 1. We will refer 
to this approximation as linear approximation in later 
discussions, and the validity of this approximation will 
become clear later. It should be noted that the nonlinear- 
ity of the coupled sine-Gordon equation is still retained in 
the equation for Pf, Eq. ([24]) . With this approximation, 
we can calculate Ai 



Ai 



ik\ — iijp- — (iuj ' 



(25) 



/(f (1 - CA^^^) exp(iP,'') cos(fcia;)da;. 



where Pi — ^ 

The IV characteristic is given by the current conser- 
vation 



(fc^-w2)2+^2^2' 



(26) 



where the first term at the right-hand side is the nor- 
mal current and the second term is the dc part of the 
Josephson current. 
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The equation for Pi is given by 

^^/'f-^cos(fc,.)A(^)expHPr). (27) 

Equation ([27|) has many solutions, such as the trivial 
vacua solution and solutions with phase kink, which will 
be discussed separately in the following subsections. 

A. State with kink 

Equation (j27p has solutions with (2m+l)7r kink with m 
being an integer. [30I.I3H Let us consider the two simplest 
periodic configurations of P[ depicted in Figs. [2la) and 
(b) where Pf ~ fiP'^^ with /; = ±1 depending on I, 
which diagonalize Eq. ([77]) 

dlP'° = CgRe(^i) cos(A:i2;) sin P"°, (28) 

where q = I for the configuration in Fig. [2ja) and q = 2 
for the configuration in Fig. [H^b). It should be noted that 
other periodic configurations are also possible. Equation 
(|28p is invariant under the transformation x L ~ x 
and P"^ <— {2m + l)7r — P'^^ , which clearly renders a kink 
at the center of junction. Equation (|28p subject to the 
boundary condition dxP'^^ = is solved numerically and 
the results are detailed in the Fig. [21 where the (2m+l)7r 
phase kink of characteristic length Ap = 1/ ■\/C(7|Re(Ai)| 
is at the center of junction x = L/2. It is this (2m + 
l)7r phase kink that pumps the dc power into plasma 
oscillation. 

In Fig. m we can see that d^P^^ forms an unquantized 
static vortex with characteristic length Ap, therefore 
it doesn't contribute to the net supercurrent (sinP;)a;t. 
There is a dc magnetic field in each layer associated with 
the static vortex. As it points in opposite directions in 




Voltage per layer [mV] 

FIG. 3: (Color online). IV characteristics calculated by 
numerical simulations (symbols) and the linearized theory 
(lines). The inset is an enlarged view. 
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FIG. 4: (Color online), (a) Radiation power at the first cur- 
rent step obtained by simulations (symbol), the linear the- 
ory (solid line) and the power balance relation Eq. (I33p 
(dashed line). The results are obtained with C — 0.000177 
and R — 707.1. The arrow is the starting point of the first 
current step. 



different junctions, the total magnetic field across the 
intrinsic Josephson junctions vanishes. Therefore it is 
impossible to realize the kink state in a single junction. 

In addition to the (2m + l)7r phase kink, there ex- 
ist the well known solitons with 2t: phase variation su- 
perposing to the (2m + l)7r phase kink, as observed in 
our simulations (not shown in Fig. [2|). In the region of 
|a;o — L/2\ » Ap, because cos(fcia;) is almost a constant 
in the narrow region A' = 1/ Cq\Ke{Ai) cos(fcia:o)|, Eq. 
(PS)) can be approximated as 

X'^dlP'° =smP'°. (29) 

Equation (P^)) has the usual soliton solution P**^ = 
4arctan[exp((x — a;o)/A')]. The total resultant Pf is 
(2m + l)7r phase kink at the center of junction and soli- 
tons with ±2tt phase variation away from the center of 
junction. The solitons with ±27r phase variation don't 
contribute to Fi as well as the net supercurrent, and 
therefore are omitted in the following discussions. 

The IV characteristics shown in Fig. [3] is calculated 
from Eq. (|26p. One remarkable feature in the /l^ charac- 
teristics is the self-induced current step, i.e., IV branch 
with constant voltage. From Eq. it is found that 

the IV characteristics for kink solutions with different 
m, e.g., the kink solutions in Fig. [5] (c),(e),(g) and (i), is 
almost the same, because the kink renders itself approxi- 
mately as a step function and the dc current contributed 
from the kink region of width Ap is negligible. It should 
be noted that there exist two branches at the cavity res- 
onance, and the right branch has a negative differential 
resistance. 

The radiation power can be readily calculated from Eq. 
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(I16p . With the linear approximation, the power is 

5r = cos6i|Aicj|V2|^|. (30) 

The dependence of on 9 andJ_Z| is consistent with the 
results shown in Fig. 5 of Ref . [30] . The results of are 
displayed in Fig. ID Similar to the IV characteristics, the 
radiation powers for the states with different phase-kink 
configurations are almost the same. 

To check the applicability of the analytical treatment, 
we solve the equation of motion Eq. (fTTj) by computer 
simulations. [301] The time step in all simulations is set to 
At = 0.0018 and the mesh size is set to Ax = 0.002. 
The accuracy is checked with smaller At and Ax. We 
use the periodic boundary condition along the z axis to 
minimize the surface effect, and attach an effective RC 
circuit to the junctions as the boundary condition along 
the X direction. |5Q|, [KJ In this case, the impedance is 
Z = R — i/Cu!, where R is the resistance and C is the 
capacitance of the RC circuit. R and C are chosen to 
make sure that \Z\ >> 1. 

The simulation results of the IV characteristics and 
radiation power are presented in Figs. O and 21 For the 
IV characteristics, there is a good agreement between 
the theory and simulation, except that the theory is in- 
capable of describing the height of the current step. For 
the radiation power, the linear theory is valid off/near 
resonance but fails inside the current steps. The fail- 
ure of the analytical treatment is caused by the strong 
plasma oscillation and existence of harmonics in the cur- 
rent step.[30t 

To derive a better estimate of the power radiation at 
the current steps, we resort to the power balance equation 
which is valid in the whole region of the IV characteris- 
tics. Taking the plasma solution Eq. ([TT]) and substitut- 
ing into the power balance equation Eq. (I14p . we obtain 
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FIG. 5: (Color online). Radiation power on current steps. 
The symbols are from the simulations and the lines are given 
by Eq. (|33j. The results are obtained with C = 0.000177 and 
R = 707.1. 



the IV characteristics in the presence of radiation 

-Jet = + f E ^^^^^f + Tri\T. (^-^^■)^- (31) 

i=i ' ' i=i 

In other words, the radiation power can be evaluated 
if we know the IV characteristics. Therefore it is useful 
to introduce an effective conductance defined as 

Je,t = ={I3 + Pd + A>, (32) 

00 

where f3d = j J2 U^j)^ conductance due to damp- 

00 

ing of plasma oscillation, /3i. = ^ U^j)^ lti6 con- 

i=i 

ductance due to radiation at both edges. From Eq. (|3ip . 
we can calculate Aj from the IV characteristics. The 
radiation power at one edge then can be evaluated by 

5, = A.u;V2 = J,c./(-^-t-|), (33) 
z cos L 

where Jo = Joxt — f3i^ is the excess current. From the fore- 
going analysis, if we can screen the radiation at one edge, 
the radiation at the other edge is enhanced. It should be 
remarked that not the whole energy pumped into plasma 
oscillation radiates into outside space. Most part of it 
is damped by dissipations inside the intrinsic Josephson 
junctions. The radiation power at current steps obtained 
by numerical simulations and power balance condition 
is depicted in Fig. [5] In contrast to the linear approx- 
imation, the estimated radiation power by Eq. p3p is 
consistent with simulations even inside the current steps 
where the amplitude of plasma oscillation is large and 
high harmonic components are present. The power in- 
creases linearly with the Joxt and the maximum power 
is as high as 8000W/cm^ from simulations (at the 6th 
cavity mode) . The maximal total radiation power at the 
first cavity mode is about lOmW if we use a mesa of sim- 
ilar dimension as the experiments [27| , which is capable 
of practical applications. 

The cavity quality factor Qc at the cavity resonance 

= fci is given by 

^ _^ Energy Stored ^ uj 
^"-^ Power Loss /3 -f 4cos6l/L|Z| ' ^ ' 

which has the order of magnitude of 100 for (3 — 0.02 
and \Z\ » 1. The half- width of the radiation frequency 
spectrum T = uj/Qc is about lOGHz, so that the radia- 
tion is almost monochromatic. The efficiency defined as 
the ratio of the radiation power to the total power input 
is 

Je m , 1 
*5e = 1-/(7— ^ + T- (35) 

Jcxt 4 cos ti L 

The efficiency Qe at the current step corresponding to a 
lower cavity mode is larger than that of a higher mode 
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FIG. 6: (Color online). Configurations of Jasym, Jsym and electromagnetic wave at successive time ti — O.OGTi, t2 — O.SlTi, 
t3 = O.SSTi and t4 = 0.94Ti, where Ti = 2ii/k\ is the period of plasma oscillation at the first cavity mode. We have subtracted 
the static part of The top figures are taken at the bottom of the first current step while the bottom figures are at the top 
of the first current step in Fig. [3] 



because of the smaller ohmic dissipations. Q^. at the top 
of the first current step in Fig. [5] is about 7.5%. 

The distribution of the c-axis uniform EM wave along 
the X direction of the junctions, as well as the supercur- 
rent, obtained from simulations at the bottom and top of 
the first current step are shown in Fig. [SI The supercur- 
rent has the same period as Pi along the c axis, and we 
only visualize it at one layer. We divide the Josephson 
current into the symmetric part Jgym and antisymmet- 
ric part Jasym with respect to the center of junction, i.e., 
sinP/(a;) = Jasym (a;) + Jsym (a;). Off resonance, Jsym is 
zero except the center of junction, while Jasym oscillates 
from —1 to +1; the magnetic field is symmetric and elec- 
tric field is antisymmetric with respect to the center of 
junction. However, at the top of the current step, the 
even part of Josephson current becomes more important, 
so the radiation is of dipole type. The corresponding 
distribution for EM wave is neither symmetric nor an- 
tisymmetric because the higher harmonics in Eq. (fT7)) 
become important. 



B. State without kink 

Equation (|77|) also has trivial vacua solutions P"" — 
mn. Without losing generality, we take m = 0. From 
Eq. (pS)) . we know that the transverse plasma cannot 
exist without a kink. Therefore the solution becomes 



Pi = Lut — iAi exp(zci;i). 



(36) 



which is nothing but the McCumber solution. Here only 
the fundamental mode is taken, which is sufficient be- 



cause of the small plasma oscillation in this solution. 
Then Eq. (gS]) is reduced to 



(37) 



and its corresponding IV characteristics without radia- 
tion is 



J b: 



l3io- 



f3 



(38) 
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FIG. 7: (Color online). Radiation power from the state with- 
out kink and its corresponding IV characteristics. Symbols 
are for simulations and lines are for theory. The vertical 
dashed line is the retrapping point obtained from the theory. 
The inset is the frequency spectrum at the strongest radia- 
tion. The results are obtained with C = 0.716 and R = 10.0. 
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FIG. 8: (Color online). Configurations of phase P, magnetic field and electric field -E^ at successive time ti = O.OTTr, 
t2 = O.lTTr, ta = 0.33Tr and t4 = OATi, where Tr is the period at the retrapping point. The phase is normalize into [0, 2n] and 
we have subtracted the static part of . The results are obtained by simulations with C — 0.716 and R = 10.0. 



The radiation power at one edge obtained with Eq. (|T6l 



Sr =cose/[2{uj^ + f3^)\Z\ 



(39) 



From the power balance condition, the IV characteristics 
with radiation is given by 



cos 6 



(40) 



where the last term represents the correction due to ra- 
diation. The minimum value of Joxt in Eq. (|40|) is the 
retrapping current Jj-, at which the input power becomes 
insufficient for the phase particle to travel across the 
damped tilted washboard potential. In the weak damp- 
ing limit /3 << 1 as in the present system, Jj- is 



3/4 



3cos6' 
1Z\L 



3/3 

Y 



1/4 



(41) 



Its corresponding voltage is oji- — (1.5-1- 
3cos9/\Z\Lpy^^ > 1, which justifies the approxi- 
mation made in Eqs. p6|) and (|37p. 

The IV characteristics calculated from the analytic 
formula and numerical simulations are presented in Fig. 
[71 A good agreement between the simulation and theory 
can be spotted. This further verifies the approximation of 
neglecting the effect of radiation on the phase dynamics 
in junctions when Z is large. The radiation power in- 
creases continuously with decreasing current and reaches 
the maximum at the retrapping point. The local minima 
in the curve are caused by the small spatial modulation 
of electromagnetic field, which changes with the voltage, 
similar to a cavity behavior. The frequency harmonics 
is undiscernible even at the maximum radiation because 
the plasma oscillation is weak. The distributions for P, 
and obtained by numerical simulations with open 
boundary condition are displayed in Fig. [51 where there is 



a small phase gradient created by radiation which is hard 
to see in the present scale. The magnetic field is antisym- 
metric with respect to the center of the junction, while 
the electric field is almost uniform along the x direction, 
except the small modulation created by radiation. 

The state without kink (McCumber state) is unstable 
in the certain region of IV curve. In a long Josephson 
junction, the system evolves into soliton states due to 
the parametric instability, [s^ We have investigated the 
stability of the state without kink in a stack of intrinsic 
Josephson junctions. The system favors the state with 
kink due to the instability of the state without kink near 
the cavity resonance. 



IV. ENERGETIC ANALYSIS 

Similar to conventional laser systems, most part of the 
input power is stored and dissipated in the junctions 
and only a small portion radiates into space. Therefore 
it is worthy of looking at the energy oscillation inside 
the junctions. One might consider that the state with 
kink costs more energy than that without kink. To see 
whether the state with kink can be realized in reality, 
it is necessary to know the energy cost to construct the 
kink. In this section, we calculate the energy stored in 
the intrinsic Josephson junctions. 

As can be read from Eq. ([7]), the system energy con- 
sists of the magnetic energy Eb , electric energy Ee and 
Josephson coupling 

Eb = (9.P^M-i5,P),t/2iV, 
Ee - {dtP^dtP).t/2N, 

= (^(l-cosPi)),t/iV, 
I 

where the energy has been normalized by the num- 
ber of layers N. In the state with kink, the mag- 
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FIG. 9: (Color online). Energy per junction stored in a thick stack of intrinsic Josephson junctions in the state with and 
without kink, (a) magnetic, (b) electric, (c) Josephson, (d) total energy. The results are obtained without radiation, and the 
results for the state with kink are obtained in the first current step shown in Fig. |3l Because the energy for different kink 
configurations is same, only the energy for the kink configuration P*"^ — [— vr, +tv, — tt, +7r] is shown in this figure. 



netic energy has contribution from the static kink 
Ebs = {d:,P''^M-'^d^P')^t/'2N and from the plasma 



oscillation = {d^P M~^dxP)xt/'2N . Here we 

show that Ebs « £^Bp- From Eq. d^Pf 
has the order of magnitude y/ (^q\Ke{Ai)\ in the nar- 
row region l/-\/C'?|R-e(^i)|. Thus, the order of E^s is 
^y\qRe{Ai)\/( << 1. On the other hand, Ebp is propor- 
tional to {Aiki)^, which is of order of 10 with the pa- 
rameters used in the present system. This also indicates 
that the magnetic energy for different kinks is roughly 
the same (thus we only show the magnetic energy for 
one kink configuration in Fig. UKa)). It is quite different 
from the usual solitons in a single Josephson junction. 

With the linear expansion of Josephson current, similar 
to Eq. ([26|l. the Josephson energy Ej in the state with 
kink can be obtained 



E, 



1 



{kl-Lo')\F,\yA 

(fc2 - cj2)2 + /32^2 



(42) 



It first decreases and then increases inside the current 
step, while is close to unity off resonance. The calculation 
of the electric energy and total energy is straightforward. 
The results are shown in Figs. [9] (b) and (d). The to- 
tal energy for different kinks is approximately the same. 
Therefore the states with kink occupy finite volumes in 
the phase space with the same energy, which makes this 
state easily accessible. 

In the state without kink, the magnetic energy is ob- 
viously zero. The Josephson coupling is 



1 



1 



2(cj2 + /32) 



(43) 



It decreases from its maximum at the retrapping point 
and saturates to unity at large currents, which is consis- 
tent with the results shown in Fig. ^c). Ej in the region 
of Jext < "^r is very close to 0, which cannot be described 
by Eq. (j43p because the system is retrapped into super- 
conducting state. The total energy is the same as that 



of state with kink in the linear ohmic IV curve, while it 
is larger than that of state with kink in the current step. 

The electric energies in Fig. [3] is much larger than 
other energy in the state without kink in a sharp con- 
trast to cases at equilibrium, where different energy con- 
tributions are expected to be the same. The reason can 
be understood if we consider Eq. pT|) as the equation 
of motion of phase particle in a titled washboard poten- 
tial. In the presence of external current, the phase par- 
ticle is accelerated and start to run in the tilted wash- 
board potential. In response to the modulated potential, 
small oscillations of phase particle are created in addition 
to the motion with a constant velocity. Meanwhile, the 
motion of phase particle causes dissipation. The steady 
state is reached when the input power and dissipation 
are balanced. On the other hand, the magnetic energy 
and Josephson coupling are solely contributed from the 
small oscillation of the phase particle. As the most part 
of input power converted into the motion with a constant 
velocity, the electric energy occupies the most part of en- 
ergy stored in the system, as shown in Fig. [5) However, 
in the state with kink, as a significant portion of the input 
power is converted into plasma oscillation in the current 
steps, rather than to solely increase the velocity of the 
phase particle, the electric energy and magnetic energy 
become comparable to each other. 



V. RADIATION PATTERN 

In this section, we calculate the far-field radiation pat- 
tern for the mesa operated in the state with kink and 
without kink, which is important both for applications 
and for differentiating various states. 

To calculate the radiation pattern, we resort to the 
Huygens principle in which the pattern is determined by 
the oscillation of the electromagnetic fields at the edges 
of samples, which can be casted into the edge magnetic 
current and electric current in the formula of equivalence 
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FIG. 10: (Color online), (a) Coordinate system for the calculation of radiation pattern from the mesa, (b) Radiation pattern 
of the mode (1, 0). (c) Radiation pattern of the mode (1, 1). Here Lx = 80^m, Ly = 300^m and Lz = 1/xm. 



principle. [53| Since there exists a significant impedance 
mismatch, the electric current produced by the oscillat- 
ing magnetic field is much smaller than the magnetic 
current produced by the oscillating electric field, so we 
can neglect the contribution from the electric current. 
The equivalent magnetic current Mg in the dimension- 
less units is given by 



(44) 



where Ee is the oscillating electric field at the edges of 
mesa and the vector n is normal to the edges. As we 
know that the radiation pattern critically depends on the 
geometry of the source, we need to consider the 3D sys- 
tem. The extension from previous analysis of 2D system 




FIG. 11: (Color online). Radiation pattern of the state with- 
out kink biased at the retrapping point. Here Lx — SO/im, 
Ly = SOOfim and Lz = l^m. The anisotropy of the pattern 
in the xy plane is due to the fact Ly » L^. 



system are sketched in Fig. [TUTa). We use the similar 
dimension as in the experiments. p?!. [2^ i.e. = 80/xm, 
Ly = 300/^m and Lz = 1/xm. Because k^^L^ << 1 with 
= oj/c, the sources at different z coordinates do not 
interfere too much, and can be treated as uniform. In 
this case, the far-field Poynting vector in the dimension- 
less units is 



cj'Li 



327r2r2e3/2 



\G\h 



(45) 



with 



G = / Afe(r')exp(-i^r' • er)(er x eiO^l', (46) 

"'edges 

where the integral is taken over the perimeter of the 
crystal. [H, [s^l With the size we use, the interference is 
mainly contributed from the source along the y direction 
since Ly is comparable to the wavelength. 

In the state with kink, the oscillation of the electric 
field in the frequency domain can be well described by 



AiLo coii{nxi: / L,j,) cos{nyTT / Ly 



(47) 



for the cavity mode (ux, riy) when the plasma oscillation 
is weak. [33| The radiation pattern from the mode {nx,ny) 
can be evaluated with Eq. (|45|) and is reported in Refs. 
[H, [13 • Inside the current step, the higher harmonics 
become important so numerical simulations are needed. 
We use computer simulations to calculate the oscillation 
of electric field at edges and then substitute the results 
into Eqs. (pij) and (PSj) to obtain the radiation pattern. 
The results at current steps corresponding to the cavity 
modes (1, 0) and (1, 1) are shown in Figs, [ini(b) and (c). 
For mode (1,0), the radiation power is maximal at the 
top of the mesa 6 = 0, and it has a maximum at the 
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FIG. 12: (Color online). Radiation power from the zero-field 
steps caused by soliton motions. The inset is the IV char- 
acteristics. The vertical dashed lines are the assignment of 
cavity mode according to fcn = rni/L with n being an integer 
and the ac Josephson relation. The results are obtained with 
C = 0.000177 and R = 707.1. 



middle of while a minimum at the middle of Ly ; for 
mode (1, 1), the radiation power is minimal at the top of 
the mesa = 0, and at the middle of and Ly. 

In the state without kink, the oscillation of the electric 
field is homogenous in the xy plane, which corresponds to 
the (0, 0) mode. The radiation pattern at the retrapping 
point is presented in Fig. 1111 It has a minimum both 
at the top of the mesa and at the middle of Ly, and a 
maximum at the middle of Lx- The anisotropy of the 



50000 



40000 



20000 



10000 - 




pattern in the xy plane is due to the fact Ly » L^ 



VI. STATE WITH SOLITONS 

To be comprehensive, we present here the results of nu- 
merical simulations on the state with solitons. It should 
be noted that in the present system, the length scale is 
Ac rather than Aj as in conventional Josephson junctions 
and in the presence of magnetic field. Therefore, to have 
solitons, the length of the junctions must be larger than 
Ac- Because of repulsive interaction, it is believed to be 
hard to achieve in-phase motion of solitons in a stack of 
junctions, despite some simulations suggest that the soli- 
tons in high velocity have attractive interaction, (ssl . [ssj 
Here we investigate the radiation due to soliton motions 
in a single junction, which is equivalent to a stack of 
junctions if one realizes the in-phase motion of solitons 
in different jimctions. 

It is well known that periodic motions and reflections 
of solitons and anti-solitons give birth to the zero-field 
steps 'At V = 2n7r/L, which corresponds to the even 
cavity modes, [s^, [57| with n the number of solitons. 
When a soliton hits the boundary, it emits an electro- 
magnetic pulse, fssi Here we only investigate the radiation 
from zero-field steps and will not discuss the Cherenkov 
radiation, [sil We perform computer simulations to trace 
out all the zero-field steps. We use L = 5Ac so that there 
exist five steps. The IV characteristics is shown in the 
inset of Fig. [T^l where the zero-field steps occur at the 
voltage corresponding to the even cavity modes. The ra- 
diation power at each step is shown in Fig. 1121 which is 
higher than that from the state with kink if one assumes 
in-phase motion of solitons. The discontinuous drops in 
the radiation power when ramping up the current are 
caused by the change in the wavelength of the Josephson 
plasma excited by the motion of solitons. The frequency 
spectrum at the first zero- field step is sketched in Fig. [13] 
and there are many frequency harmonics with the fun- 
damental frequency not satisfying the ac Josephson rela- 
tion. In the state with solitons, the radiated frequency 
depends on the configuration of solitons except the first 
zero-field step.jg^l It is noticed that, as indicated in Fig. 
[T51 the fundamental frequency and voltage at all the steps 
never satisfy the ac Josephson relation, contrasting with 
the state of kink revealed theoretically 30] and the exper- 
imental observations 1271 [2911 . 
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VII. DISCUSSIONS AND CONCLUSIONS 



FIG. 13: (Color online). Frequency spectrum at the first 
zero-field step at Jcxt = 0.57. The vertical dashed line is 
the frequency given by the ac Josephson relation with voltage 
V — 2n/L. The fundamental peak does not come to the 
dashed line means that the ac Josephson relation is broken. 
The results are obtained with C = 0.000177 and R = 707.1. 



In the present work, the phase dynamics and its elec- 
trodynamics in a thick stack of intrinsic Josephson junc- 
tions in the absence of external magnetic field are in- 
vestigated both analytically and numerically. There is 
good consistency between the analytical theory and sim- 
ulations. 
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In the state with phase kink, there are many current 
steps at both even and odd cavity modes. The phase kink 
plays a role of coupUng the plasma to the cavity modes, 
as such the plasma oscillation is largely enhanced. The 
radiation power from the state with phase kink is 
times larger than that without phase kink. The plasma 
oscillation is uniform through the c axis. Thus the far- 
field radiation power grows as N squared, the so-called 
superradiation. At the bottom of the first current step, 
the magnetic field is symmetric with respect to the center 
of the junction. The antisymmetric component becomes 
more and more important when going into the current 
step. The states with different phase kink configurations 
are degenerate in the sense that they have the same IV 
characteristics, power radiation and energy stored in the 
system. 

In the state without phase kink, the plasma oscilla- 
tion docs not couple to the cavity mode, and thus the 
radiation power is very small. The power increases with 
decreasing current and reaches the maximum at the re- 
trapping point. The radiation occurs in a broad region 
of voltage. The frequency satisfies the ac Josephson rela- 
tion, and high frequency harmonics are almost invisible. 
The magnetic field is antisymmetric in this state. The 
far-field radiation pattern in this state is quite different 
from that in the state with phase kink, which is a clear 
fingerprint of the dynamic state realized by the system. 

In the state with solitons, electromagnetic pulses are 
radiated from junctions when solitons hit the boundary. 
There are many frequency harmonics, but the fundamen- 
tal frequency never satisfies the ac Josephson relation. It 
would be ideal for exciting strong terahertz wave with 
solitons because the power is about 25000 W/cm^, pre- 
suming one could realize the in-phasc motion of solitons 
in a thick stack of long intrinsic Josephson junctions. 

It is illuminating to discuss the dynamic state realized 
in the recent experiments for terahertz radiation [13, 
in light of the present theoretical analysis. Coherent 
radiations were detected in the resistive curve in Ref. 
[27| with the frequency corresponding to the first cavity 
mode. One order of the magnitude stronger radiations 
were observed in the region of voltage with anomalous IV 
characteristics in Ref. [29| , and there are many frequency 
harmonics at a given voltage. In both experiments, the 
frequency obeys the ac Josephson effect and thus the 
state with traveling solitons can be ruled out. On the 
other hand, large cavity resonances cannot be excited 
in the state without phase kink; furthermore, the radi- 
ation from the state without phase kink occurs weakly 
in a wide range of voltage. Therefore, it is unlikely rel- 



evant to the experimental observations. The state with 
phase kink, in contrast, seems to be consistent with the 
experiments so far. In this state, the plasma oscillation 
is uniform through the stack of Josephson junctions, it 
thus supports superradiation as observed in the experi- 
ments. Moreover, the periodic arrangement of static kink 
along the stack direction allows to pump dc powers into 
large plasma oscillations, which yields self-induced cur- 
rent steps. It is noticed that, to obtain the overall shape 
of the IV curve, we need to take the heating effect into 
account. More works are needed to clarify the synchro- 
nization process. 

It should be remarked that there are propagating waves 
besides the standing wave in Eq. P7|) , because of the ra- 
diation. The amplitude of the propagating waves has 
the order of magnitude of l/\Z\ and thus can be safely 
neglected as the first-step approximation, which is con- 
firmed by the numerical simulations. The radiation has 
only some negligible effects on the dynamics inside the 
junctions, which permits us to calculate the radiation 
perturbatively. When the mismatch of impedance at the 
edges is reduced, e.g., the thickness of a stack of intrin- 
sic Josephson junctions is comparable to Ac, one has to 
consider the radiation and interference for the analysis of 
phase dynamics inside the junctions self-consistently. 

The essential property of the stack of junctions for 
realization of the state with kink is the strong induc- 
tive coupling. As for other possible effects on the state 
with phase kink, we find in the simulations that this 
state is very stable against small magnetic fields, thermal 
fluctuations. [6l| In Ref. [3l|, it is shown that the state 
with kink is stable against the modulation of critical cur- 
rent. Thus it is likely to be realized experimentally. The 
state with kink is promising for the application of tera- 
hertz radiation. It is also useful for terahertz detectors, 
amplifiers and mixers. 
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